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Abstract 

Protein translation is a multistep process which can be represented as a cascade of biochemi- 
cal reactions (initiation, ribosome assembly, elongation, etc.), the rate of which can be regulated by 
small non-coding microRNAs through multiple mechanisms. It remains unclear what mechanisms of 
microRNA action are most dominant: moreover, many experimental reports deliver controversal mes- 
sages on what is the concrete mechanism actually observed in the experiment. Parker and Nissan [37j 
demonstrated that it is impossible to distinguish alternative biological hypotheses using the steady 
state data on the rate of protein synthesis. For their analysis they used two simple kinetic models 
of protein translation. In contrary, we show that dynamical data allow to discriminate some of the 
mechanisms of microRNA action. We demonstrate this using the same models as in |37j for the sake of 
comparison but the methods developed (asymptotology of biochemical networks) can be used for other 
models. As one of the results of our analysis, we formulate a hypothesis that the effect of microRNA 
action is measurable and observable only if it affects the dominant system (generalization of the limit- 
ing step notion for complex networks) of the protein translation machinery. The dominant system can 
vary in different experimental conditions that can partially explain the existing controversy of some of 
the experimental data. 



1 Introduction 

MicroRNAs (miRNAs) arc currently considered as key regulators of a wide variety of biological path- 
ways, including development, differentiation and oncogenesis. Recently, remarkable progress was made 
in the understanding of microRNA biogenesis, function and mechanism of action. Mature microRNAs 
are incorporated into the RISC effector complex, which includes as a key component an Argonaute 
protein. MicroRNAs affect gene expression by guiding the RISC complex toward specific target mR- 
NAs. The exact mechanism of this inhibition is still a matter of debate. In the past few years, several 
mechanisms have been reported, some of which are contradictory (for review, see [5J O EI])- These 
include in particular inhibition of translation initiation (acting at the level of cap-40S or 40S-AUG- 
60S association steps), inhibition of translation elongation or premature termination of translation. 
MicroRNA-mediated mRNA decay and sequestration of target mRNAs in P-bodies have been also 
proposed. Moreover, some microRNAs mediate target mRNA cleavage [56], chromatin reorganization 
followed by transcriptional repression or activation [26j [41] , or translational activation [40] [50] . 
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The most frequently reported, but also much debated, mechanism of gene repression by microRNAs 
occurs at the level of mRNA translation. This is particularly true when the microRNA is not fully 
complementary to its target. At this level, several mode of actions have been suggested (see Fig. [T]). 
Historically, the first proposed mechanism was inhibition of translation elongation. The major argu- 
ment supporting this hypothesis was the observation that the inhibited mRNA remained associated 
with the polysomal fraction (in which mRNAs are associated with polysomes), an observation that 
was reproduced in different systems [TSJ [351 E21 13S1 ISS] ■ The idea of a post-initiation mechanism was 
further supported by the observation that some mRNAs can be repressed by a microRNA even when 
their translation is cap-independent (mRNAs with an IRES or A-capped) [lj [23 E3J [43] . Although 
it was initially proposed that the ribosomes was somehow "frozen" on the mRNA, it is important 
to note that it is difficult to discriminate experimentally between different post-initiation potential 
mechanisms, such as elongation inhibition, premature ribosome dissociation ("ribosome drop-off') or 
normal elongation but nascent polypeptide degradation. The last proposition (which can occur in 
conjonction with the two others) is supported by the fact that the mRNA-polysomal association is 
puromycin-sensitive, indicating that it depends on a peptidyl-transferase activity [2| I35j. However, no 
nascent peptide has ever been experimentally demonstrated; thus degradation would occur extremely 
rapidly after synthesis [l2"l 135] . Degradation, in any case, is proteasome-independent [32] . Prema- 
ture ribosome dissociation is supported by decreased read-through of inhibited mRNA [43] . Ribosome 
drop-off and/or ribosomal "slowing" are supported by the slight decrease in the number of associated 
ribosomes observed in some studies J38J El] . 

Concurrently, several reports have been published indicating an action of microRNAs at the level 
of initiation. An increasing number of papers reports that microARN-targetted mRNAs shift toward 
light fractions in polysomal profiles [31 [23 S3 • This shows a decrease of mature translating ribosomes, 
suggesting that microRNAs act on the initiation step. Moreover, several reports show that microRNA- 
mediated inhibition is relieved when translation is driven by a cap-independent mechanism such as 
IRES-mRNA or A-capped- mRNA [2Ql [27l [42]. This observation was confirmed in several in- vitro 
studies 36, 48, 53, 51.. In particular, in one of those, an excess of eIF4F could relieve the inhibition, 
and inhibition lead to decreased 80S in polysomal gradient [55] . 

Most of the data indicating a shift toward light polysomal fraction or requirement for a cap- 
dependent translation are often interpreted in favour of involvement of microRNAs at early steps of 
translation, i.e., cap binding and 40S recruitment. However, some of them are also compatible with 
a block at the level of 60S subunit joining. This hypothesis is also supported by in-vitro experiments 
showing a lower amount of 60S relative to 40S on inhibited mRNAs. Moreover, toe-printing experiments 
show that 40S is positioned on the AUG [52]. Independently, it was shown that eIF6, an inhibitor of 
60S joining, is required for microRNA action [7], but this was contradicted by other studies [9]. 

Thus, data regarding the exact step of translational inhibition arc clearly contradictory. Taking 
also into account the data about mRNA degradation and P-bodies localization, it is difficult to draw a 
clear picture of the situation, and the exact mechanism by which microRNA repress mRNA expression 
is highly controversial, let alone the interrelations between the different mechanisms and their possible 
concomitant action. Several attempts to integrate the different hypothesis have been made [6l [9l ITTT [2TT 
[29J [49] . For example, one mechanism could act as a "primary" effect, and the other as a "secondary" 
mechanism, either used to reinforce the inhibition or as back-up mechanism. In others, the different 
mechanisms could all coexist, but occur differentially depending on some yet unidentified characteristic. 
For example, it has been observed than the same mRNA targeted by the same microRNA can be 
regulated either at the initiation or the elongation step depending on the mRNA promoter and thus 
on mRNA nuclear history [28] . It was also proposed that technical (experimental) problems, including 
the variety of experimental systems used, may also account for these discrepancies [6l El [UJ . However, 
this possibility does not seem to be sufficient to provide a simple and convincing explanation to the 
reported discrepancies. 

A possible solution to exploit the experimental observations and to provide a rational and straight- 
forward interpretation is the use of mathematical models for microRNA action. Nissan and Parker 
recently studied the steady states of two simple kinetic models |37j and provided critical analysis for 
the experiments with alternative mRNA cap structures and IRES elements [36l [48l [51] , leading to 
possible explanation of the conflicting results. The authors suggested that the relief of translational 
repression upon replacement of the cap structure can be explained if microRNA is acting on a step 
which is not rate-limiting in the modified system, in which case, an effect of microRNA can simply 
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Figure 1: Interaction of microRNA with protein translation process. Four mechanisms of translation 
repression which are considered in the mathematical modeling are indicated: 1) on the initiation process, 
preventing assembling of the initiation complex; 2) on a late initiation step, such as searching for the 
start codon; 3) on the ribosome assembly; 4) on the translation process. There exist other mechanisms 
of microRNA action on protein translation (transcriptional, transport to P-bodies, ribosome drop-off, 
co-translational protein degradation and others) that are not considered in this paper. Here 40S and 
60S are light and heavy components of the ribosome, 80S is the assembled ribosome bound to mRNA, 
eIF4F is an translation initiation factor, PABC1 is the Poly-A binding protein, "cap" is the mRNA cap 
structure needed for mRNA circularization, RISC is the RNA-induced silencing complex. 

not be observed. It was claimed that it is impossible to discriminate between two alternative interpre- 
tations of the biological experiments with cap structure replacement from the sole monitoring of the 
steady state level of protein [37] . 

Two remarks can be made in this regard. First, in practice not only the steady state level of 
protein can be observed but also other dynamical characteristics, such as the relaxation time, i.e. the 
time needed to achieve the steady state rate after a perturbation (such as restarting the translation 
process). We argue that having these measurements in hands, one can distinguish between two alterna- 
tive interpretations. In this paper we provide such a method from the same models as constructed by 
Nissan and Parker, for comparison purposes. However, the method applied can be easily generalized 
for other models. 

Second, even in the simple non-linear model of protein translation, taking into account the recycling 
of ribosomal components, it is difficult to define what is the rate limiting step. It is known from the 
theory of asymptotology of biochemical networks [17] that even in complex linear systems the "rate 
limiting place" notion is not trivial and can not be reduced to a single reaction step. Moreover, in 
non- linear systems the "rate limiting place" can change with time and depend on the initial conditions. 
Hence, conclusions of 37] should be re-considered for the non-linear model, made more precise and 
general. The notion of rate limiting step should be replaced by the notion of dominant system. 

In this paper we perform careful analysis of the Nissan and Parker's models and provide their 
approximate analytical solutions, which allows us to generalize the conclusions of [37] and make new 
testable predictions on the identifiability of active mechanism of microRNA-dependent protein trans- 
lation inhibition. 
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2 Methods 



2.1 Asymptotology and dynamical limitation theory for biochemical reac- 
tion networks 

Most of mathematical models that really work are simplifications of the basic theoretical models and 
use in the backgrounds an assumption that some terms are big, and some other terms are small enough 
to neglect or almost neglect them. The closer consideration shows that such a simple separation on 
"small" and "big" terms should be used with precautions, and special culture was developed. The 
name "asymptotology" for this direction of science was proposed by [30] defined as "the art of handling 
applied mathematical systems in limiting cases" . 

In chemical kinetics three fundamental ideas were developed for model simplification: quasiequilib- 
rium asymptotic (QE), quasi steady-state asymptotic (QSS) and the idea of limiting step. 

In the IUPAC Compendium of Chemical Terminology (2007) one can find a definition of limiting 
step [44^ : "A rate-controlling (rate-determining or rate- limiting) step in a reaction occurring by a 
composite reaction sequence is an elementary reaction the rate constant for which exerts a strong 
effect - stronger than that of any other rate constant - on the overall rate." 

Usually when people are talking about limiting step they expect significantly more: there exists 
a rate constant which exerts such a strong effect on the overall rate that the effect of all other rate 
constants together is significantly smaller. For the IUPAC Compendium definition a rate-controlling 
step always exists, because among the control functions generically exists the biggest one. On the 
contrary, for the notion of limiting step that is used in practice, there exists a difference between 
systems with limiting step and systems without limiting step. 

During XX century, the concept of the limiting step was revised several times. First simple idea 
of a "narrow place" (the least conductive step) could be applied without adaptation only to a simple 
cycle or a chain of irreversible steps that are of the first order (see Chap. 16 of the book [23] or the 
paper by [5] or the section |2~51 of this paper). When researchers try to apply this idea in more general 
situations they meet various difficulties. 

Recently, we proposed a new theory of dynamic and static limitation in multiscale reaction networks 
|16[ I17j . This approach allows to find the simplest network which can substitute a multiscale reaction 
network such that the dynamics of the complex network can be approximated by the simpler one. 
Following the asymptotology terminology [54j . we call this simple network the dominant system (DS). 
In the simplest cases, the dominant system is a subsystem of the original model. However, in the 
general case, it also includes new reactions with kinetic rates expressed through the parameters of the 
original model, and rates of some reactions are renormalized: hence, in general, the dominant system 
is not a subsystem of the original model. 

The dominant systems can be used for direct computation of steady states and relaxation dynamics, 
especially when kinetic information is incomplete, for design of experiments and mining of experimental 
data, and could serve as a robust first approximation in perturbation theory or for preconditioning. 
Dominant systems serve as correct generalization of the limiting step notion in the case of complex 
multiscale networks. They can be used to answer an important question: given a network model, 
which are its critical parameters? Many of the parameters of the initial model are no longer present 
in the dominant system: these parameters are non-critical. Parameters of dominant systems indicate 
putative targets to change the behavior of the large network. 

Most of reaction networks are nonlinear, it is nevertheless useful to have an efficient algorithm for 
solving linear problems. First, nonlinear systems often include linear subsystems, containing reactions 
that are (pseudo)monomolecular with respect to species internal to the subsystem (at most one in- 
ternal species is reactant and at most one is product). Second, for binary reactions A + B — > if 
concentrations of species A and B (ca, cb) are well separated, say > cb then we can consider this 
reaction as B — > ... with rate constant proportional to ca which is practically constant, because its 
relative changes are small in comparison to relative changes of cb ■ We can assume that this condition 
is satisfied for all but a small fraction of genuinely non-linear reactions (the set of non-linear reactions 
changes in time but remains small). Under such an assumption, non- linear behavior can be approxi- 
mated as a sequence of such systems, followed one each other in a sequence of "phase transitions" . In 
these transitions, the order relation between some of species concentrations changes. Some applica- 
tions of this approach to systems biology are presented by [45] . The idea of controllable linearization 



4 



"by excess" of some reagents is in the background of the efficient experimental technique of Temporal 
Analysis of Products (TAP), which allows to decipher detailed mechanisms of catalytic reactions [551 . 

Below we give some details on the approaches used in this paper to analyse the models by Nissan 
and Parker [37]. 

2.2 Notations 

To define a chemical reaction network, we have to introduce: 

• a list of components (species); 

• a list of elementary reactions; 

• a kinetic law of elementary reactions. 

The list of components is just a list of symbols (labels) A\, ...A n . Each elementary reaction is repre- 
sented by its stoichiometric equation 

si si 

where s enumerates the elementary reactions, and the non- negative integers a 8 i, j3 s i are the stoichio- 
metric coefficients. A stoichiomentric vector j s with coordinates "/ S i = /3 s j — a S i is associated with 
each elementary reaction. 

A non-negative real extensive variable iVj > 0, amount of Ai, is associated with each component 
Ai. It measures "the number of particles of that species" (in particles, or in moles). The concentration 
of Ai is an intensive variable: Cj = Ni/V, where V is volume. In this paper we consider the volume 
(of cytoplasm) to be constant. Then the kinetic equations have the following form 

^i = J2w s (c,T) ls + v, (2) 

s 

where T is the temperature, w s is the rate of the reaction s, v is the vector of external fluxes normalized 
to unite volume. It may be useful to represent external fluxes as elementary reactions by introduction 
of new component together with income and outgoing reactions — » Ai and Ai — > 0. 

The most popular kinetic law of elementary reactions is the mass action law for perfect systems: 

W s (c,T)=k s [[cf" i , (3) 

where k s is a "kinetic constant" of the reaction s. 



2.3 Quasy steady-state and quasiequilibrium asymptotics 

Quasiequilibrium approximation uses the assumption that a group of reactions is much faster then 
other and goes fast to its equilibrium. We use below superscripts ' f ' and ' s ' to distinguish fast and 
slow reactions. A small parameter appears in the following form 

^ = £ <(c,TK + - £ ]T ^(c,T) 7? f , (4) 

s, slow ^, fast 

To separate variables, we have to study the spaces of linear conservation law of the initial system ([4]) 
and of the fast subsystem 

fast 

If they coincide, then the fast subsystem just dominates, and there is no fast-slow separation for 
variables (all variables are either fast or constant). But if there exist additional linearly indepen- 
dent linear conservation laws for the fast system, then let us introduce new variables: linear func- 
tions 6 1 (c), ...b n (c), where 6 1 (c), ...b m (c) is the basis of the linear conservation laws for the initial 
system, and & 1 (c), ...b m+l (c) is the basis of the linear conservation laws for the fast subsystem. Then 
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b m+l+ (c), ...b n {c) are fast variables, b m+1 (c), ...b m+l (c) are slow variables, and fo 1 (c), ...b m (c) are con- 
stant. The quasiequilibrium manifold is given by the equations ^ ? w\(c, T)j* = and for small e it 
serves as an approximation to a slow manifold. 

The quasi steady-state (or pseudo steady state) assumption was invented in chemistry for descrip- 
tion of systems with radicals or catalysts. In the most usual version the species are split in two groups 
with concentration vectors c s ( "slow" or basic components) and c ( "fast intermediates" ) . For catalytic 
reactions there is additional balance for c f , amount of catalyst, usually it is just a sum bf — J2i c i- The 
amount of the fast intermediates is assumed much smaller than the amount of the basic components, 
but the reaction rates are of the same order, or even the same (both intermediates and slow components 
participate in the same reactions). This is the source of a small parameter in the system. Let us scale 
the concentrations c f and c s to the compatible amounts. After that, the fast and slow time appear and 
we could write c s = W s (c s , c f ), c f = (c s , c f ), where e is small parameter, and functions W s , W l are 
bounded and have bounded derivatives (are "of the same order" ) . We can apply the standard singular 
perturbation techniques. If dynamics of fast components under given values of slow concentrations is 
stable, then the slow attractive manifold exists, and its zero approximation is given by the system of 
equations W i {c s 1 c i ) = 0. 

The QE approximation is also extremely popular and useful. It has simpler dynamical proper- 
ties (respects thermodynamics, for example, and gives no critical effects in fast subsystems of closed 
systems). Nevertheless, neither radicals in combustion, nor intermediates in catalytic kinetics are, 
in general, close to quasiequilibrium. They are just present in much smaller amount, and when this 
amount grows, then the QSS approximation fails. 

The simplest demonstration of these two approximation gives the simple reaction: S + E <-> SE — * 
P + E with reaction rate constants kf 1 and ki . The only possible quasiequilibrium appears when the 
first equilibrium is fast: kf = n^/e. The corresponding slow variable is C s — cs+cse, bE — ce+cse = 

const. For the QE manifold we get a quadratic equation j^-cse — cs c e — (C s ~ c SE)(bE ~ cse)- 

1 

This equation gives the explicit dependence cse(C s ), and the slow equation reads C s = —k2CsE(C s ), 
C s + cp — bs = const. 

For the QSS approximation of this reaction kinetics, under assumption 6e <C bs, we have fast 
intermediates E and SE. For the QSS manifold there is a linear equation k^ cgCE — k^ csE~k 2 csE = 0; 
which gives us the explicit expression for cse(cs)'- c se — &i Cs&e/(^i~ c S + k^ + A; 2 ) (the standard 
Michaelis-Menten formula). The slow kinetics reads cs — —k^cs(bE — cse(cs)) + c Se{cs)- The 
difference between the QSS and the QE in this example is obvious. 

The terminology is not rigorous, and often QSS is used for all singular perturbed systems, and QE is 
applied only for the thermodynamic exclusion of fast variables by the maximum entropy (or minimum 
of free energy, or extremum of another relevant thermodynamic function) principle (MaxEnt). This 
terminological convention may be convenient. Nevertheless, without any relation to terminology, the 
difference between these two types of introduction of a small parameter is huge. There exists plenty 
of generalizations of these approaches, which aim to construct a slow and (almost) invariant manifold, 
and to approximate fast motion as well. The following references can give a first impression about these 
methods: Method of Invariant Manifolds (MIM) ([HUB], Method of Invariant Grids (MIG), a discrete 
analogue of invariant manifolds ([13]), Computational Singular Perturbations (CSP) ([3TJ [351 [57] ) 
Intrinsic Low-Dimensional Manifolds (ILDM) by [33], developed further in series of works by [3]), 
methods based on the Lyapunov auxiliary theorem (|24j). 

2.4 Multiscale monomolecular reaction networks 

A monomolecular reaction is such that it has at most one reactant and at most one product. Let us 
consider a general network of monomolecular reactions. This network is represented as a directed graph 
(digraph) [UJ: vertices correspond to components A i} edges correspond to reactions A4 — > Aj with 
kinetic constants kji > 0. For each vertex, Ai, a positive real variable Cj (concentration) is defined. 
"Pseudo-species" (labeled 0) can be defined to collect all degraded products, and degradation reactions 
can be written as Ai — > with constants fcoj. Production reactions can be represented as — > Ai 
with rates kio. The kinetic equation for the system is 

— kio + ^ t kijCj — ^ ^ kjiCi, (5) 
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or in vector form: c = Kq + Kc. Solution of this system can be reduced to a linear algebra problem: 
let us find all left and right (r % ) eigenvectors of K , i.e.: 

Kr i = X i r i ,l i K = X i l i , (6) 
with the normalisation r l ) = Sij, where <Jy is Kronecker's delta. Then the solution of (|5|) is 

n 

c(t) = c s + r k (l k ,c(0) - c s ) exp(-A fc t), (7) 
fe=l 

where c s is the steady state of the system ((SJ), i.e. when all ^ = 0, and c(0) is the initial condition. 

If all reaction constants fey would be known with precision then the eigenvalues and the eigenvectors 
of the kinetic matrix can be easily calculated by standard numerical techniques. Furthermore, Singular 
Value Decomposition (SVD) can be used for model reduction. But in systems biology models often 
one has only approximate or relative values of the constants (information on which constant is bigger 
or smaller than another one). Let us consider the simplest case: when all kinetic constants are very 
different (separated), i.e. for any two different pairs of indices I = J = {i',f) we have either 

k] ^> kj or kj *C fcj. In this case we say that the system is hierarchical with timescales (inverses of 
constants fey, j ^ 0) totally separated. 

Linear network with totally separated constants can be represented as a digraph and a set of orders 
(integer numbers) associated to each arc (reaction). The lower the order, the more rapid is the reaction. 
It happens that in this case the special structure of the matrix K (originated from a reaction graph) 
allows us to exploit the strong relation between the dynamics ([5]) and the topological properties of the 
digraph. In this case, possible values of U k are 0, 1 and the possible values of r% are -1, 0, 1 with 
high precision. In previous works, we provided an algorithm for finding non-zero components of h, Ti, 
based on the network topology and the constants ordering, which gives a good approximation to the 
problem solution [TB1[T71H5] . 

2.5 Dominant system for a simple irreversible catalytic cycle with limiting 
step 

A linear chain of reactions, A\ — > — > ...An, with reaction rate constants ki (for Ai — > A;+i), 
gives the first example of limiting steps. Let the reaction rate constant k q be the smallest one. Then 
we expect the following behaviour of the reaction chain in time scale > l/k q : all the components 
A%, ...A q -i transform fast into A q , and all the components A q +i, ...A n -\ transform fast into A n , only 
two components, A q and A n are present (concentrations of other components are small) , and the 
whole dynamics in this time scale can be represented by a single reaction A q A n with reaction 
rate constant k q . This picture becomes more exact when k q becomes smaller with respect to other 
constants. 

The kinetic equation for the linear chain is 

c% — ki~iCi—i kiCi, (8) 

The coefficient matrix K of these equations is very simple. It has nonzero elements only on the main 
diagonal, and one position below. The eigenvalues of K are — ki (i — 1, ...n — 1) and 0. The left and 
right eigenvectors for eigenvalue, 1° and r°, are: 

Z° = (l,l,...l), r° = (0,0,...0,l), (9) 

all coordinates of l a are equal to 1, the only nonzero coordinate of r° is r° and we represent vector- 
column r° in row. 

The catalytic cycle is one of the most important substructures that we study in reaction networks. 
In the reduced form the catalytic cycle is a set of linear reactions: 

Ai A 2 -> ...A n 

Reduced form means that in reality some of these reaction are not monomolecular and include 
some other components (not from the list A\, . . . A n ). But in the study of the isolated cycle dynamics, 
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concentrations of these components are taken as constant and are included into kinetic constants of 
the cycle linear reactions. 

For the constant of elementary reaction A; — > we use the simplified notation ki because the product 
of this elementary reaction is known, it is A+i for i < n and A\ for i = n. The elementary reaction 
rate is Wi — kid, where Ci is the concentration of A{. The kinetic equation is: 

Ci = fcj_iCj_! - fcjCj, (10) 

where by definition cq = c n , ko = k n , and wq — w n . In the stationary state (ci = 0), all the Wi are 
equal: Wi — w. This common rate w we call the cycle stationary rate, and 

w = -j— - -; Ci = — , (11) 

IT + • ' ■ fcl 1 

where b — Ci is the conserved quantity for reactions in constant volume. Let one of the constants, 
fcmin, be much smaller than others (let it be £; m in = k n ): 



h > fctnin if i^n . (12) 



In this case, in linear approximation 



k 

d = b-^-, w = k n b. (13) 



The simplest zero order approximation for the steady state gives 

c n = b, cj=0(t^n). (14) 

This is trivial: all the concentration is collected at the starting point of the "narrow place" , but may 
be useful as an origin point for various approximation procedures. 

So, the stationary rate of a cycle is determined by the smallest constant, £; m j n , if it is much smaller 
than the constants of all other reactions (fT2"|) : 

w w k min b. (15) 

In that case we say that the cycle has a limiting step with constant k m { n . 

There is significant difference between the examples of limiting steps for the chain of reactions and 
for irreversible cycle. For the chain, the steady state does not depend on nonzero rate constants. It 
is just c„ = b, c\ — C2 = ... = c„_i = 0. The smallest rate constant k q gives the smallest positive 
eigenvalue, the relaxation time is r = l/k q . The corresponding approximation of eigenmode (right 
eigenvector) r 1 has coordinates: r\ = ... = r] l _ 1 = 0, r x q — 1, r* +1 = ... = r\_ x = 0, r n = — 1. 
This exactly corresponds to the statement that the whole dynamics in the time scale > l/k q can be 
represented by a single reaction A q — ► A n with reaction rate constant k q . The left eigenvector for 
eigenvalue k q has approximation I 1 with coordinates l\ = l\ = ... = l q = 1, l q+1 = ... = Z„ = 0. 
This vector provides the almost exact lumping on time scale > l/k q . Let us introduce a new variable 
ciump = J2i li c i> i- e - c iump = c\ + C2 + •■■ + Cq- For the time scale > l/k q we can write ci ump + c n » 6, 
dci um p/dt ~ ^qCi um p, dc n / dt ~ k q c\ um p. 

In the example of a cycle, we approximate the steady state, that is, the right eigenvector r° for 
zero eigenvalue (the left eigenvector is known and corresponds to the main linear balance b: q = 1). 
In the zero-order approximation, this eigenvector has coordinates = ... = r„_i = 0, rjj = 1. 

If k n jki is small for all i < n, then the kinetic behaviour of the cycle is determined by a linear chain 
of n— 1 reactions A\ — > A^ — > ...A n , which we obtain after cutting the limiting step. The characteristic 
equation for an irreversible cycle, n"=i(A + — Tii=i K i = tends to the characteristic equation for 
the linear chain, AJ^™^ (A + ki) = 0, when k n — > 0. 

The characteristic equation for a cycle with limiting step (knjki <C 1) has one simple zero eigenvalue 
that corresponds to the conservation law ^2 a = b and n — 1 nonzero eigenvalues 

Xi = -h + Si{i< n). (16) 
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where Si — > when Y,i< n i~ ~~ * °- 

A cycle with limiting step (|10[) has real eigenspectrum and demonstrates monotonic relaxation 
without damped oscillations. Of course, without limitation such oscillations could exist, for example, 
when all fci = fc > 0, (i = 1, ...n). 

The relaxation time of a stable linear system (|10|) is, by definition, r = 1/ min{i?e(— A;)} (A ^ 0). 
For small k n , r w l/fc T , &t = rninjfci}, (i = 1, ...n — 1). In other words, for a cycle with limiting step, 
k T is the second slowest rate constant: fc m ; n «C k T < .... 

3 Results 

3.1 Model assumptions 

We consider two models of action of microRNA on protein translation process proposed in [37] : the sim- 
plest linear model, and the non-linear model which explicitly takes into account recycling of ribosomal 
subunits and initiation factors. 

Both models, of course, represent significant simplifications of biological reality. First, only a lim- 
ited subset of all possible mechanisms of microRNA action on the translation process is considered (see 
Fig.[T]). Second, all processes of synthesis and degradation of mRNA and protein are deliberately ne- 
glected. Third, interaction of microRNA and mRNA is simplified: it is supposed that when microRNA 
is added to the experimental system then only mRNA with bound microRNAs are present (this also 
assumes that the concentration of microRNA is abundant with respect to mRNA). Concentrations of 
microRNA and mRNA are supposed to be constant. Interaction of only one type of microRNA and 
one type of mRNA is considered (not a mix of several microRNAs). The process of initiation is greatly 
simplified: all initiation factors are represented by only one molecule which is marked as eIF4F. 

Finally, the classical chemical kinetics approach is applied, based on solutions of ordinary differential 
equations, which supposes sufficient and well-stirred amount of both microRNAs and mRNAs. Another 
assumption in the modeling is the mass action law assumed for the reaction kinetic rates. 

It is important to underline the interpretation of certain chemical species considered in the system. 
The ribosomal subunits and the initiation factors in the model exist in free and bound forms, moreover, 
ribosomal subunits can be bound to several regions of mRNA (the initiation site, the start codon, the 
coding part). Importantly, several copies of fully assembled ribosome can be bound to one mRNA. 
To model this situation, we have to introduce the following quantification rule for chemical species: 
amount of "ribosome bound to mRNA" means the total number of ribosomes translating proteins, 
which is not equal to the number of mRNAs with ribosome sitting on them, since one mRNA can hold 
several translating ribosomes (polyribosome). In this view, mRNAs act as places or catalyzers, where 
translation takes place, whereas mRNA itself formally is not consumed in the process of translation, 
but, of course, can be degraded or synthesized (which is, however, not considered in the models 
described further). 

3.2 The simplest linear protein translation model 

The simplest representation of the translation process has the form of a circular cascade of reactions 
[37] (see Fig. |2). 

The list of chemical species in the model is the following: 

1. 40S, free small ribosomal subunit. 

2. mRNA:40S, small ribosomal subunit bound to the iniation site. 

3. AUG, small ribosomal subunit bound to the start codon. 
The catalytic cycle is formed by the following reactions: 

1. 40S — > mRNA:40S, Initiation complex assembly (rate fcl). 

2. mRNA:40S — > AUG, Some late and cap-independent initiation steps, such as scanning the 
5'UTR by for the start codon AUG recognition (rate k2). 

3. AUG — > 40S, combined processes of 60S ribosomal unit joining and protein elongation, which 
leads to production of the protein (rate fc3), and fall off of the ribosome from mRNA. 

The model is described by the following system of equations [37] : 
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( d [405] (t) 
dt 



= -fci[40S] + k 3 [AUG] 



d [mRNA : 405] (t) 
Jt 

d[AUG](t) 



= fci [405] - fc 2 [mi?AA : 405] 



(17) 



dt 



= k 2 [mRNA : 405] - k 3 [AUG] 



, Prsynth{t) = k 3 [AUG](t) 

where Prsynth(t) is the rate of the protein synthesis. 

Following [57] . let us assume that £3 3> fci, fc 2 . This choice was justified by the following statement: 
"...The subunit joining and protein production rate (k3) is faster than kl and k2 since mRNA40S com- 
plexes bound to the AUG without the 60S subunit are generally not observed in translation initiation 
unless this step is stalled by experimental methods, and elongation is generally thought to not be rate 
limiting in protein synthesis.." |37j . 

Under this condition, the equations (|17[) have the following approximate solution (which becomes 
the more exact the smaller the kl k k2 ratio): 



[405] (t) 

[mRNA : 405] (i) 
[AUG](t) 



for the initial condition 



[405]q 

— + — 

fci ^ fc 2 



" 1/fcl " 


1 

+ k 3 


' -1 " 


l/*2 


1 


l/*3 








-k 3 t 



1 
~2 





1 

-1 



-(fci+fca)* 



Prsynth(t) — -j 














' [405] 




"[405] " 




[mi?AA : 405] 









[AUG] 









Prsynth 





0. 



(18) 
(19) 



(20) 



From the solution (fig)) it follows that the dynamics of the system evolves on two time scales: 1) fast 
elongation dynamics on the time scale « ^/k 3 ; and 2) relatively slow translation initiation dynamics 
with the relaxation time t re i ss k ± +k 2 ' protein synthesis rate formula (|19p does not include the 
k 3 rate, since it is neglected with respect to ki,k 2 values. From (| 19[) we can extract the formula for 
the protein synthesis steady-state rate Prsynth (multiplier before the parentheses) and the relaxation 
time trei for it (inverse of the exponent power): 



Prsynth 



[405] 



(21) 



Now let us consider two experimental situations: 1) the rates of the two translation initiation steps 
are comparable k\ w fc 2 ; 2) the cap-dependent rate k\ is limiting: k\ <C fc 2 - Accordingly to [37], 
the second situation can correspond to modified mRNA with an alternative cap-structure, which is 
much less efficient for the assembly of the initiation factors, 40S ribosomal subunit and polyA binding 
proteins. 

For these two experimental systems (let us call them 'wild-type' and 'modified' correspondingly), 
let us study the effect of microRNA action. We will model the microRNA action as a mechanism which 
effectively diminish the value of a kinetic rate coefficient. Let us assume that there are two alternative 
mechanisms: 1) microRNA acts in a cap-dependent manner (thus, reducing the k\ constant) and 
2) microRNA acts in a cap-indcpcndcnt manner, for example, through interfering with 60S subunit 
joining (thus, reducing the /c 2 constant). The dependence of the steady rate of protein synthesis 
Prsynth ~ j_ and the relaxation time t re i as kl \. k2 on the efficiency of the microRNA action (i.e., 

how much it is capable to diminish a rate coefficient) is shown on Fig. [3l 

Interestingly, experiments with cap structure replacement were made and the effect of microRNA 
action on the translation was measured J36[ 08] . No change in the protein rate synthesis after applying 
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Figure 2: The simplest model of microRNA action on the protein translation, represented with use 
Systems Biology Graphical Notation (a) and schematically with the condition on the constants (b). T 
two mechanism of microRNA action. 



microRNA was observed. From this it was concluded that microRNA in this system should act through 
cap-dependent mechanism (i.e., the normal 'wild- type' cap is required for microRNA recruitement). It 
was argued that this could be a misinterpretation [37] since in the 'modified' system, cap-dependent 
translation initiation is a rate limiting process (k\ «C k^), hence, even if microRNA acts in the cap- 
independent manner (inhibiting A^), it will have no effect on the final steady state protein synthesis 
rate. They confirmed this by the graph similar to the Fig. [3ji. 

From the analytical solution (fl8|) we can further develop this idea and claim that it is still possible to 
detect the action of microRNA in the 'modified' system if one measures the protein synthesis relaxation 
time: if it significantly increases then microRNA probably acts in the cap-independent manner despite 
the fact that the steady state rate of the protein synthesis does not change (see the Fig. [3b ) . This is 
a simple consequence of the fact that the relaxation time in a cycle of biochemical reactions is limited 
by the second slowest reaction [T7J . If the relaxation time is not changed in the presence of microRNA 
then we can conclude that none of the two alternative mechanisms of microRNA-based translation 
repression is activated in the system, hence, microRNA action is dependent on the structure of the 
'wild-type' transcript cap. 

The observations from the Fig. [3J are recapitulated in the Table [TJ This analysis (of course, over- 
simplified in many aspects) provides us with an important lesson: observed dynamical features of 
the translation process with and without presence of microRNA can give clues on the mechanisms of 
microRNA action and help to distinguish them in a particular experimental situation. Theoretical 
analysis of the translation dynamics highlights what are the important characteristics of the dynamics 
which should be measured in order to infer the possible microRNA mechanism. 

This conclusion suggests the notion of a kinetic signature of microRNA action mechanism which 
we define as the set of measurable characteritics of the translational machinery dynamics (features 
of time series for protein, mRNA, ribosomal subunits concentrations) and the predicted tendencies of 
their changes as a response to microRNA action through a particular biochemical mechanism. 



3.3 The non-linear protein translation model 

To explain the effect of microRNA interference with translation initiation factors, a non-linear version 
of the translation model was proposed [37 which explicitly takes into account recycling of initiation 
factors (eIF4F) and ribosomal subunits (40S and 60S). 

The model contains the following list of chemical species (see also Fig. |4j: 

1. 40S, free 40S ribosomal submit. 

2. 60S, free 60S ribosomal subunit. 

3. eIF4F, free initiation factor. 
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Figure 3: Graphs illustrating the predicted change in the steady-state rate of protein synthesis (left), 
and its relaxation time (i.e., the time needed to recover from a perturbation to the steady state value), 
right. Four curves are presented. The black ones are for the wild-type cap structure, which is modeled by 
k% = &2- The red ones are for the modified structure, when fci < The main conclusion from the left 
graph is that if microRNA acts on a late initiation step, diminishing k<i then it's effect is not measurable 
unless &2 is very strongly suppressed (as reported in |37j). The main conclusion from the right is that 
the effect of microRNA should be measurable in this case if one looks at dynamical features such the 
relaxation time. 



Table 1: Modeling of two mechanisms of microRNA action in the simplest linear model. microRNA 
action effect is described for the protein synthesis steady rate and the relaxation time. It is assumed that 
the ribosome assembly+elongation step in protein translation, described by the k^ rate constant, is not 
rate limiting 



Observable value 


Initiation (ki) 


Step after initiation, cap-independent (^2) 


Elongation (£3) 


Steady-state rate 
Relaxation time 

Steady-state rate 
Relaxation time 


decreases 
increases slightly 

decreases 
no change 


Wild-type cap 

decreases 
increases slightly 
A- cap 

no change 
increases drastically 


no change 
no change 

no change 
no change 
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4. mRNA:40S, formed initiation complex (containing 40S and the initiation factors), bound to the 
initiation site of mRNA. 

5. AUG, initiation complex bound to the start codon of mRNA. 

6. 80S, fully assembled ribosome translating protein. 

There are four reactions in the model, all considered to be irreversible: 

1. 40S + eIF4F — > mRNA:40S, assembly of the initiation complex (rate fci). 

2. mRNA:40S — -> AUG, some late and cap-independent initiation steps, such as scanning the 5'UTR 
by for the start codon AUG recognition (rate fo). 

3. AUG — > 80S, assembly of ribosomes and protein translation (rate £3). 

4. 80S — > 60S+40S, recycling of ribosomal subunits (rate fc 4 ). 
The model is described by the following system of equations [37] : 



d [405] 
dt 

d [elFAF 
dt 

d \mRNA : 405 



= -fc x [405] [elFAF] + fc 4 [805] 



= -fei[405] [elFAF] + k 2 [mRNA : 405] 



dt 

d [AUG] 

~dt " 
d [605] 



dt 
d [805] 



= fci[405] [elFAF] - k 2 [mRNA : 405] 
k 2 [mRNA : 405] - k 3 [AUG] [605] 
k 3 [AUG] [605] + fc 4 [805] 



(22) 



k 3 [AUG] [605] - fc 4 [805] 
k 3 [AUG] [605] 



dt 

^ Prsynth(t) 

where [405] and [605] are the concentrations of free 40S and 60S ribosomal subunits, [elFAF] is a 
concentration of free translation initiation factors, [mRNA : 405] is the concentration of 40S subunit 
bound to the initiation site of mRNA, [AU G] is the concentration of the initiation complex bound to 
the start codon, [805] is the concentration of ribosomes translating protein, and Prsynth is the rate 
of protein synthesis. 

The model (|22l) contains three independent conservations laws: 



iRNA : 405] + [405] + [AUG] + [805] = [405] , 



(23) 



[mRNA : 405] + [elFAF] = [eIF4F] , 

[605] + [805] = [605] o, 
The following assumptions on the model parameters were suggested in |37] : 



(24) 
(25) 



fc 4 <C fci, ha, k 3 - k 3 > h, k 2 - [e/F4F] <C [405] ; [e/F4F] < [605] < [405] (26) 

with the following justification: "...The amount 40S ribosomal subunit was set arbitrarily high ... as it 
is thought to generally not be a limiting factor for translation initiation. In contrast, the level of eIF4F, 
as the canonical limiting factor, was set significantly lower so translation would be dependent on its 
concentration as observed experimentally... Finally, the amount of subunit joining factors for the 60S 
large ribosomal subunit were estimated to be more abundant than eIF4F but still substoichiometric 
when compared to 40S levels, consistent with in vivo levels... The k4 rate is relatively slower than the 
other rates in the model; nevertheless, the simulations overall protein production was not altered by 
changes of several orders of magnitude around its value... " |37j . 

Notice that further in our paper we show that the last statement about the value of fc 4 is needed 
to be made more precise: in the model by Nissan and Parker, fc 4 is a critical parameter. It does not 
affect the steady state protein synthesis rate only in one of the possible scenarios (inefficient initiation, 
deficit of the initiation factors). 
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3.3.1 Steady state solution 



The final steady state of the system can be calculated from the conservation laws and the balance 
equations among all the reaction fluxes: 



k 2 ■ [mRNA : 405] s = k 3 ■ [AUG] S ■ [605] s = k 4 ■ [805] s = k x ■ [405] s • [eIF4F]„, (27) 
where 's' index stands for the steady state value. Let us designate a fraction of the free [60S] ribosomal 

[60S] S 
[60S]o 1 



subunit in the steady state as x = [Sri j° . Then we have 



[mRNA : 405] s = ^[605] (1 - a;), [AUG] S = ^ f— 
k 2 k 3 \ x 

[eIF4F] s = [eIF4F} - £[60S] (1 - x), 

k 2 (28) 

[605] s = [605] a;, [805] s = [605] (1 - x), 
[405], = [405] - [605] (1 - x)(l + ^) - ^ ( ' ' 



k 2 k 3 \ x 

and the equation to determine x, in which we have neglected the terms of smaller order of magnitude, 
based on conditions 



x 3 +x 2 {a + (S - 1) + 09 - 1)) +x (-a + (6 - l)(/3 - 1)) + 7(1 -0)=O, 

k 2 [eIF4F] k 2 k 4 [405] ■ (29) 



fci[605] [605] fc 4 ' ' /c 3 [605]o' [605] 

From the inequalities on the parameters of the model, we have S > 1, 7 <C 1 and, if ki 3> 
k4/[eIF4F]o then a <C /3. From these remarks it follows that the constant term 7(1 — (3) of the 
equation (|29|) should be much smaller than the other polynomial coefficients, and the equation (|29|) 
should have one solution close to zero and two others: 



Xq 



fc4 

&3-([405]o-[605] )' 

„ aP [e!FAF] k 2 kj[eIFAF] 1 

Xl P+ S~p [605] fc 4 + fc 1 fc 4 [405] 1 - K F nTt 2 ' (30) 



2^2 



[405] o 
[605] : 



provided that a<C|l — 5|ora<C|l — (3\. In the expression for x\ we can not neglect the term 
proportional to a, to avoid zero values in (|29|) . 

The solution x 2 is always negative, which means that one can have one positive solution xq <C 1 if 
^[dot^aki 2 w ^ an< ^ tw0 positive solutions xo an d xi if ~ 1- However, from {28]), (|30|) and 

(1261) it is easy to check that if x\ > then so does not correspond to positive value of [eIF4F] s . This 
means that for a given combination of parameters satisfying (|26p we can have only one steady state 
(either xq or x\). 

The two values x = xq and ir = x\ correspond to two different modes of translation. When, for 
example, the amount of the initiation factors [eIFAF]o is not sufficient to provide efficient initiation 
(\eIFAF]o < [gogf fc 4 ) x = x i) then most of the 405 and [605] subunits remain in the free form, 
the initiation factor [elFAF] being always the limiting factor. If the initiation is efficient enough 
{[eIFAF]o > [ 6n g| nfc4 )i then we have x = xq <C 1 when almost all 605 ribosomal subunits are engaged 
in the protein elongation, and [elFAF] being a limiting factor at the early stage, however, is liberated 
after and ribosomal subunits recycling becomes limiting in the initiation (see the next section for the 
analysis of the dynamics). 

Let us notice that the steady state protein synthesis rate under these assumptions is 
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miRNA 




Figure 4: Non-linear model of microRNA action on the protein translation, represented with use of 
Systems Biology Graphical Notation (a) and schematically with the condition on the constants (b). 
The difference from the simplest model (|2|) is in the explicit description of initiation factors eIF4F, and 
ribosomal subunits 40S and 60S recycling. 



f, [RnQ1 ., [eIF4F] k 2 
Prsynth = k± ■ [606] (1 - x) » < [606] fc 4 . (31) 

[ fc 2 ■ [eIFAF] , else 

This explains the numerical results obtained in |37j : with low concentrations of [eIF4F] microRNA 
action would be efficient only if it affects fc 2 or if it competes with elFAF for binding to the mRNA cap 
structure (thus, effectively further reducing the level [eIF4F]o). With higher concentrations [eIF4F]o, 
other limiting factors become dominant: [606]o (availability of the heavy ribosomal subunit) and k± 
(speed of ribosomal subunits recycling which is the slowest reaction rate in the system). Interestingly, 
in any situation the protein translation rate does not depend on the value of k\ directly (of course, 
unless it does not become "globally" rate limiting), but only through competing with elFAF (which 
makes the difference with the simplest linear protein translation model). 

Equation pip explains also some experimental results reported in |36j : increasing the concentration 
of [elFAF] translation initiation factor enhances protein synthesis but its effect is abruptly saturated 
above a certain level. 

It would be interesting to make some conclusions on the shift of polysomal profile from the steady 
state solutions [30] In this model, the number of ribosomes sitting on mRNA N po i ysome is defined by 
^polysome = [mHiVA] ' wnere [mRNA] is the concentration of mRNA. However, [mRNA] is not an 
explicit dynamical variable in the model, it is implicitly included in other model constants, such as 
ki, together with the effective volume of cytoplasmic space considered in the model. Nevertheless, the 
model can predict the relative shift of the polysome profile. In the steady state 

K4 ■ 606 Li it r „„„, , > 1 

N polysome ~ 1 - x = I 41 J °' [606] fc 4 . (32) 

I k 2 ■ [eIF4F] , else 

and N po iy Some changes in the same way as the protein synthesis steady state value. 
3.3.2 Analysis of the dynamics 

It was proposed to use the following model parameters in |37| : k\ = fc 2 = 2, k% = 5, fc 4 = 1, 
[406] = 100, [606] = 25, [eIF4F]o — 6. Although all estimations made in the previous section 
remain valid for them, in this section we have modified the parameter values (anyway very arbitrary 
estimated) to make the dynamical features of the model more illustrative. After we have verified that 
the semi- analytical approach remains valid for the set of "normal" parameter values from |37j (and, 
in general, for a very large set of parameters). First, we modified fc 2 = 3 to avoid possible artifact 
symmetries related to equal parameter values. Second, we have changed fc 4 = 0.1 to have the case of 
x = xq steady state (which, as it will be shown, corresponds to the three-stage dynamics, which allows 
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to demonstrate three dominant dynamical systems instead of only one in the case of x = x\ when 
elAFA is a limiting factor during the whole relaxation process) and ks = 50 to have better separation 
between k^ and k\, &2. 

Simulations of the protein translation model with these parameters and the initial conditions 

" [40S] 

[elFAF] 

[mRNA : AOS] 

[AUG] 

[80S] 
. [60S] 

are shown on the Fig. [6] The system shows non-trivial relaxation process which takes place in several 
epochs. One can check that the "efficient translation initiation" [x = Xq <C 1) scenario is realized in 
this case. Qualitatively we can distinguish the following stages: 

1) Stage 1: Relatively fast relaxation with conditions [405] > [elFAF] , [60S] > [AUG]. During 
this stage, the two non-linear reactions 40S + elFAF -> mRNA : 40S and AUG + 60S -> 80S can 
be considered as pseudo-monomolecular ones: elFAF — * mRNA : 40S and AUG — > 80S with rate 
constants dependent on [40S] and [60S] respectively. This stage is characterized by rapidly establishing 
quasiequilibrium of three first reactions (Rl, R2 and R3 with k\, and k^ constants). Biologically, 
this stage corresponds to assembling of the translation initiation machinery, scanning for the start 
codon and assembly of the first full ribosome at the start codon position. 

2) Transition between Stage 1 and Stage 2. 

3) Stage 2: Relaxation with the conditions [40S] > [elFAF], [60S] < [AUG]. During this stage, 
the reactions 40S + elFAF -> mRNA : 40S and AUG + 60S -> 80S can be considered as pseudo- 
monomolecular elFAF — > mRNA : 40S and 60S — > 80S. This stage is characterized by two local 
quasi-steady states established in the two network reaction cycles (formed from R1-R2 and R3-R4 
reactions). Biologically, this stage corresponds to the first round of elongation, when first ribosomes 
moves along the coding region of mRNA. The small ribosomal subunit 40S is still in excess which keeps 
the initiation stage (reaction R1-R2 fluxes) relatively fast. 

4) Transition between Stage 2 and Stage 3. 

5) Stage 3: Relaxation with the conditions [40S] < [elFAF], [60S] < [AUG]. During this stage, 
the reactions 40S + elFAF -> mRNA : 40S and AUG + 60S -> 80S can be considered as pseudo- 
monomolecular 40S — > mRNA : 40S and 60S — > 80S. During this stage all reaction fluxes are 
balanced. Biologically, this stage corresponds to the stable production of the protein with constant 
recycling of the ribosomal subunits. Most of ribosomal subunits 40S are involved in protein elongation, 
so the initiation process should wait until the end of elongation for that they would be recycled. 

Stages 1-3 can be associated with the corresponding dominant systems [T7] which are shown on 
Fig. El 

Below we give a more detailed analysis of stages 1-3 and transitions between them. 

3.3.3 Stage 1: translation initiation and assembly of the first ribosome at the 
start codon 

The dominant system of the Stage 1 (Fig. [5^) can be modeled as a linear system of equations (notice 
that it is not equivalent to the system of equations that would correspond to fully monomolecular 
reaction network because the reaction R2 is still bimolecular despite the fact that the products of this 
reaction do not interact, which leads to the linear description): 



]40SJ 
[eIFAF] 



[60Slo 



(33) 
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kl'=kl 






mRNA:40S 




AUG 




80S 




k2 









kl"=kl-eIF4F k3"=k3-AUl 




Stage 1 



Stage 2 



Stage 3 



Figure 5: Dominant systems for three stages of relaxation of the model (|22j) . Stage 1) The dominant 
system is pseudo-linear network of reactions. Stage 2) The dominant system is a quasi-steady state 
approximation, where one supposes that the fluxes in two network cycles are balanced. Stage 3) The 
dominant system is pseudo- linear network of reactions. 



' d [eIF4F] 
dt 

d \mRNA : 405 



dt 



= —k[[eIF4F] + k 2 [mRNA : 405] 

= k[[eIF4:F} - k 2 [mRNA : 405] 



d [AUG] 

dt 
d [805] 

dt~ ~ 



(34) 



= k 2 [mRNA : 405] - k' 3 [AUG] 
■-k 3 [AUG] [605] - fc 4 [805] 



where k[ = ki ■ [405] , k' 3 = k 3 ■ [605] and we consider that at this stage the changes of 40S and 60S are 
relatively slow. This system has simple approximate solution, taking into account constraints on the 
parameters k 2 <C k'i,k' 3 ; fc 4 <C k[, k' 3l k 2 , also assuming k 2 <C [k' x — k' 3 \, and for the initial condition 



(35) 



[elFAF] 




" [e/F4F] - 


[mRNA : 405] 







[AUG] 







[805] 





0. 



we have 



\eIF4F](t) 
[mRNA : 405] (t) 
[AUG](t) 
.[805] 



k 2 ■ [eIFAF} 



l/*2 
l/k' 3 
l/k 4 . 



' 


1 

- 1 



-kU 



1 

^2 



-k[t ^_ 

k4 



(36) 



-k^t 



J 



From this solution, one can conclude that the relaxation of this model goes at several time scales 
(very rapid ~ min(l/k' 3 , l/k[) and slow ~ 1/^4) and that when eIF4F, mRNA:40S and AUG already 
reached their quasiequilibrium values, [80S] continues to grow. This corresponds to the quasiequilib- 
rium approximation asymptotic (see section 12. 3| . At some point 80S will reach such a value that it 
would be not possible to consider 60S constant: otherwise the conservation law ([23)) will be violated. 
This will happen when [805] << [605] condition is not satisfied anymore, i.e., following our convention 



to consider "much smaller" as difference in one order of magnitude, at t' 



[605] 



The same 



W-[eIF 4F] -fc 2 ' 

consideration is applicable for another conservation law (|23|) in which [805] is included, but from the 

[40S] _ , 
10- [elFAFj^ ± 



time point t" 



From [405] » [e/F4F] and [405] > [605] we have min(t',t") = t' . 
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Non-linear Parker's model 



Non-linear Parker's model 




time 
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g 

I io° 
o 
c 
o 






10 ' 

time 



o 

CO 
1 

I 10° 

o 

o 
o 

10 



d) 




kl*40S*eIF^ : I 
k2*mRNA: 4 0$ 
k3*AUG*60S 
k4*80S 



4 6 

time 



Figure 6: Simulation of the non-linear protein translation model with parameters k\ = 2, &2 = 3, 
&3 = 50, ^4 = 0.1, [40<S]o = 100, [605] o = 25, [e/F4F]o = 6. a) and b) chemical species concentrations at 
logarithmic and linear scales; c) and d) reaction fluxes at logarithmic and linear scales. By the dashed 
line several stages are delimited during which the dynamics can be considered as (pseudo-) linear. To 
determine where "S>" and "<C" conditions are violated, we arbitrarily consider "much bigger" or "much 
smaller" as difference in one order of magnitude (by factor 10). 
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This means that the parameters [405], [605] of the (local) steady states for [eIF4F[ and [AL/G] should 
slowly (at the same rate as [805]) change from the time point t' (variable [mRNA : 405] does not 
change because its local steady state does not depend on [405], [605]). In other words, after t = t' the 
Stage 1 solution ((36]) should be prolonged as 

[805](i) = fc2 ' l eIF/L F} ^ _ e _ fe4t ^ ^ j 40 £^ = [405]q _ [ 805 ]( t ) ! [605](i) = [605] - [805] (t), 
[eIFAF]{t) = , , [mRNA : 405](t) = [eIF4F] , [AUG](t] h ' [eIF4F]o 



fei([405] -[805](t))' L ' JW 1 JU ' L JW fc 3 ([605] -[805](i)) 

(37) 

From these equations, one can determine the effective duration of the Stage 1: by definition, it will 
be finished when one of the two conditions ([405] > [elFAF], [605] > [AUG]) will be violated, which 
happens at times ~ fc 2 .[e/F4F] an d ~ fc 2 -[e/f 4f ] corres P on dingly, hence, the second condition will be 
violated first (from [605] < [405] ). 

3.3.4 Stage 2: first stage of protein elongation, initiation is still rapid 

The Stage 2 is characterized by conditions [elFAF] < [405] , [605] < [AUG] . This fact can be used 
for deriving the quasi-steady state approximation: we assume that the reaction fluxes in two network 
cycles (R1-R2 and R3-R4) are independently balanced: 

fci[405] • [elFAF] = k 2 [mRNA : 405], k 3 [AUG] ■ [605] = fc 4 [805]. (38) 
Then (|22p is simplified and, using the conservation laws, we have a single equation on [405]: 

d[4QS](t) _ -k 2 [eIF4F] -[A0S](t) , k 4 [60S] (A-[40S](t)) 

dt ^ + [405](t) + tA +A -[40S](t) ' ( M > 

where A = [405] (i) + L4J7G] (i) is a constant quantity conserved accordingly to the quasi-steady state 
approximation. Equation (|39[) can be already integrated but let us further simplify it for our analysis. 
Having in mind fc 4 <C k% and assuming that at the beginning of the Stage 2 [>1J7G] ^> || , we can 
simplify to 

d[40S](t) _ -fc 2 [eJF4F]o-[405](t) ujanq] , AM 

dt ~ ^+[4os](t) + fc4 ^ Ulb J°' (40) 

and, further, assuming that at the beginning of the Stage 2 we have [405] ^> [j4Z7G] let us approximate 
the right-hand side of the equation by a piecewise-linear function 

d[40S](t) ^ f -k 2 [eIF4F] Q ^^^- + fc 4 [605] , if [40S](t) > ^ 
dt | _ kl [eiF4F] _ [405]^) + fc 4 [605] , [405](t) < fe, 

where [405]|t = t» is the amount of 405 at the beginning of the Stage 2. Then the descent of [405] (i) 
can be separated into linear and exponential phases: 

, ,,, ([40S]\t=t"-Ki-(t-t"),ift<t" + T 2 , x 

405 ](t) w < (421 
1 JW \ [405] s2 - ([405] s2 - [405]| t=t „) exp (-K 2 (t - t" - r 2 )) , if t > t" + r 2 ' 

where Ki,K 2 are linear and exponential slopes and [405] s2 is the quasi-steady state value of [405] at 
the end of the Stage 2: 
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t~ + [405]*=*// 
ifi = k 2 eIF4F ^ - fc 4 [605] 



K, 



2 



' 3 ^ + [405]|* 
fci [e/F4F] 



fcj 1 (43) 

t 405 ^ = ^ fca[eJ F4F] - 
fe 4 [60S]o 

[405] | t 



T2 = Ki 

Other dynamic variables are expressed through [405] (t) as 



[eIF4F](t) = [eIF4F} 



k2 

fci 



If + [405] W 

405] (t) 



[AE/G1W = [405)1,=,,, + ^ m J^^ m — ) - [mm- (44) 



[605] (t) = [605] o T 



A-3 



[80S] (i) = [605] 



§ + [AC/G](i) 
L4C/G](f) 



M + [A[/G](i) 



At some point, the amount of free small ribosomal subunit 40S, which is abundant at the beginning 
of the Stage 2, will not be sufficient to support rapid translation initiation. Then the initiation factor 
elFAF will not be the limiting factor in the initiation and the condition [405] 3> [elFAF] will be 



violated. We can estimate this time as t'' 



[40S]o 



eIF4F] ' 



3.3.5 Stage 3: steady protein elongation, speed of initiation equals to speed of 
elongation 

During the Stage 3 all fluxes in the network become balanced and the translation arrives at the steady 
state. From Fig. [5] it is clear that the relaxation goes independently in the cycle i?3 — i?4, where the 
relaxation equations are simply 



[60S](t) = [605] s + ([605]|t = t"'-[605] s )e- (fc4+fe3[ ' 4c/GI|t = t "' )(t - t "' ) ,[805](i) = [605] - [605] (t). (45) 

where t'" is the time when the Stage 3 of the relaxation starts. This relaxation goes relatively fast, 
since k^[AUG\\t=t'" is relatively big. So, during the Stage 3, one can consider the cycle Rl — i?4 
equilibrated, with [805] = [80] s , [605] = [60] s values. 

Hence, the relaxation during the Stage 3 consists in redistributing concentrations of 40S and 
mRNA:40S to their steady states in a linear chain of reactions R1 — R2 (the value of [AZ7G] is relatively 
big and can be adjusted from the conservation law (|23|) ). Using the pseudo-linear approximation of 
this stage (see Fig[5]), we can easily write down the corresponding approximate relaxation equations: 

[mRNA : 405] (t) = 

[mRNA : 40],(1 - e~ k2t ) + B( e -Me/f4F]| t=f ,,, (*-*'") _ e -fc 2 (t-t"')) + [ mRNA . 405]| t=t /»e- fe2(t -*"' ) , 
[405] (t) = [405] s + Be-^-*"") 

(46) 
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where B = (1 - ki [ eI p 4F y„ )([40] s - [405] | t=f /» )■ [ 40 ]« and ["i-RATA : 405] s are the steady-state 
values of the corresponding variables, see (|2"5)) . The values [605]|t=t'", [e/F4F]| t=t ///, L4[/G]| t=t »< and 
[miL/VA : 40S']|t = t"' can be estimated from (|4"4|) . using the [405]|t = f" value. The relaxation time at 
this stage equals 

T3 ^A^e/i^F]^ ' k 2 ' fc 4 + fc 3 [4C/G]| t=t - ' 
The solution for the Stage 3 can be further simplified if fc 2 <C fci [elFAF] \t=t"' or fc 2 3> fci [el-F^F] |t = f" . 

3.3.6 Transitions between stages 

Along the trajectory of the dynamical system (|22|) there are three dominant system each one trans- 
forming into an other. At the transition between stages, two neighbor dominant systems are united and 
then split. Theoretically, there might be situations when the system can stay in these transition zones 
for long periods of time, even infinitely. However, in the model (1221) this is not the case: the trajectory 
rapidly passes through the transition stages and jumps into the next dominant system approximation. 

Three dominant approximations can be glued, using the concentration values at the times of the 
switching of dominant approximation as initial values for the next stage. Note that the Stage 2 has 
essentially one degree of freedom since it can be approximated by a single equation (|39[) . Hence, one 
should only know one initial value [405]|t=t" to glue the Stages 1 and 2. The same is applied to the 
gluing of Stages 2 and 3, since in the end of Stage 2 all variable values are determined by the value of 

[4as]|t=f«. 

3.3.7 Case of always limiting initiation 

As it follows from our analysis, the most critical parameter of the non-linear protein translation model 
is the ratio (5 = ^j^sfo ' Above we have considered the case /3 > 1 which is characterized by a 
switch of the limiting factor in the initiation (from elFAF at the Stages 1 and 2 to 405 at the Stage 
3). 

In the case P < 1 the dynamics becomes simpler and consists of one single stage: relaxation 
accordingly to (|36|) and further with correction (|37|) with the relaxation time ~ (the quasiequilibrium 
approximation corresponding to the Stage 1 works well for the whole translation process). The reason 
for this is that if the initiation is not efficient then the system is never in the situation of the Stage 
2 conditions when the cycle R1-R2 is balanced with much bigger flux than the cycle R3-R4. This 
approximation is the more exact the smaller (3 value, however, the value of (3 should not be necessary 
very small. For example, for the default parameter values of the model (3 — 0.48, and it well reproduces 
the dynamics (see Fig. [7j>d). From numerical experiments one can see that even for = 0.95 the 
dynamics is qualitatively well reproduced. 

To model the A-cap structure effect with very weak capacity for initiation (assembly of the initiation 
factors and 40S subunit), we should also consider the case 

ki < k 4 < fc 2 < fc 3 , (47) 

for which the solution derived above is not directly applicable. However, the analytical calculations 
in this case can be performed in the same fashion as above. The detailed derivation of the solution is 
given in the Supplementary materials. The effect of putting ki very small on the steady state protein 
synthesis and the relaxation time is shown on Fig. [7l 

In a similar way all possible solutions of the equations (|22|) with very strong inhibitory effect of 
microRNA on a particular translation step can be derived. These solutions will describe the situation 
when the effect of microRNA is so strong that it changes the dominant system (limiting place of the 
network) by violating the initial constraints (|26p on the parameters (for example, by making k^ smaller 
than other kts). Such possibility exists, however, it can require too strong (non-physiological) effect of 
microRNA-dependent translation inhibition. 

3.3.8 Effect of microRNA on the translation dynamics 

As a result of the above analysis, we can assemble an approximate solution of the non-linear system 
under assumptions ([26]) about the parameters. An example of the approximate solution is given on 
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Figure 7: Examples of the exact (circles) and approximate (solid lines) solutions of the non-linear protein 
translation model, a) For our set of parameters k\ = 2, ki = 3, k% = 50, k^ = 0.1; b) For parameters 
k\ = 1, &2 = 5, /%3 = 50, /%4 = 0.01; c) For the set of parameters from [37j : d) Reaction fluxes for 
the set of parameters c). Dashed black vertical lines denote evaluated transition points between the 
dynamics stages. Dashed red vertical points denote the time points where [405"] (t) = 10 • [eIFAF]{t) and 
[405] {t) = [eIF4F](t)/W respectively. 
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Fig. [7J The advantage of such a semi-analytical solution is that one can predict the effect of changing 
the system parameters. For example, on Fig. [JJd the solution is compared to an exact numerical one, 
where the parameters have been changed but still obey the initial constraints (f26|) . 

One of the obvious predictions is that the dynamics of the system is not sensitive to variations of 
ka, so if microRNA acts on the translation stage controlled by k^ then no microRNA effect could be 
observed looking at the system dynamics (being the fastest one, k% is not a critical parameter in any 
scenario). 

If microRNA acts on the translation stage controlled by k<t (for example, by ribosome stalling 
mechanism) then we should consider two cases of efficient (j3 > 1) and inefficient (j3 < 1) initiation. In 
the first case the steady state protein synthesis rate is controlled by k± (as the slowest, limiting step) 
and any effect on k± would lead to the proportional change in the steady state of protein production. 
By contrast, in the case of inefficient initiation, the steady state protein synthesis is not affected by k^. 
Instead, the relaxation time is affected, being ~ However, diminishing k\ increases the (3 parameter, 
hence, this changes "inefficient initiation" scenario for the opposite, hence, making k± critical for the 
steady state protein synthesis anyway when k<± becomes smaller than ^Tgg~gj~~ ■ For example, for the 
default parameters of the model, decreasing k^ value first leads to no change in the steady state rate 
of protein synthesis but increases the relaxation time and, second, after the threshold value k2 ^^P° 
starts to affect the steady state protein synthesis rate directly (see Fig. UK)- This is in contradiction 
to the message from |37j that the change in k^ by several orders of magnitude does not change the 
steady state rate of protein synthesis. 

Analogously, decreasing the value of ki can convert the "efficient" initiation scenario into the 
opposite after the threshold value ^jj^fe ■ We can recapitulate the effect of decreasing &2 in the 
following way. 1) in the case of efficient initiation k^ does not affect the steady state protein synthesis 
rate up to the threshold value after which it affects it in a proportional way. The relaxation time 
drastically increases, because decreasing k^ leads to elongation of all dynamical stages durations (for 
example, we have estimated the time of the end of the dynamical Stage 2 as t'" — k2 .^p ^4F ] )• 
However, after the threshold value the relaxation time decreases together with &2, quickly dropping to 
its unperturbed value. 2) in the case case of "inefficient" initiation the steady state protein synthesis 
rate depends proportionally on the value of k^ (f3"Tj) . while the relaxation time is not affected (see 
Fig.©. 

MicroRNA action on k\ directly does not produce any strong effect neither on the relaxation time 
nor on the steady state protein synthesis rate. This is why in the original work [37] cap-dependent 
mechanism of microRNA action was taken into account through effective change of the [eIFAF] value 
(total concentration of the translation initiation factors), which is a critical parameter of the model 
(seeEl). 

The effect of microRNA on various mechanism and in various experimental settings (excess or 
deficit of eIF4F, normal cap or A-cap) is recapitulated in Table [2] The conclusion that can be made 
from this table is that all four mechanism shows clearly different pattern of behaviour in various 
experimental settings. From the simulations one can make a conclusion that it is still not possible 
to distinguish between the situation when microRNA does not have any effect on protein translation 
and the situation when it acts on the step which is neither rate limiting nor 'second rate limiting' 
in any experimental setting (fc 3 in our case). Nevertheless, if any change in the steady-state protein 
synthesis or the relaxation time is observed, theoretically, it will be possible to specify the mechanism 
responsible for it. 



4 Discussion 

The role of microRNA in gene expression regulation is discovered and confirmed since ten years, how- 
ever, there is still a lot of controversial results regarding the role of concrete mechanisms of microRNA- 
mediated protein synthesis respression. Some authors argue that it is possible that the different modes 
of microRNA action reflect different interpretations and experimental approaches, but the possibility 
that microRNAs do indeed silence gene expression via multiple mechanisms also exists. Finally, mi- 
croRNAs might silence gene expression by a common and unique mechanism; and the multiple modes 
of action represent secondary effects of this primary event [6, 9 , 11. 

The main reason for accepting a possible experimental bias could be the studies in vitro, where 
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Figure 8: Effect of decreasing some model parameters mimicking different mechanisms of miRNA action 
on translation. Relaxation time here is defined as the latest time at which any chemical species in the 
model differs from its final steady state by 10%. A) and B) correspond to the scenario with "inefficient" 
initiation, with use of the model parameters proposed in [37] {k\ = ki = 2, k% = 5, = 1, [e/-F4F]o = 6, 
^OS^o = 25, ^OS'Jo = 100), which gives (3 = 0.48 < 1. C) and D) correpond to the scenario with 
"efficient" initiation, with our choice of parameters (k\ = 2, k% = 3, k% = 50, k^ = 0.1, [e/-F4.F]o = 6, 
^OS^n = 25, ^OSJo = 100), which gives (3 = 7.2 > 1. The absciss value indicates the degree of inhibition 
(decreasing) of a parameter. E-H) same as A-D) but for a modified cap structure, modeled by reduced 
k\ parameter: k\ = 0.01 for these curves, the other parameters are the same as on A-D) correspondingly. 
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Table 2: Modeling of four mechanisms of microRNA action in the non-linear protein translation model. 
MicroRNA action effect is described for the protein synthesis steady rate and the relaxation time (see 
also Fig. [8|). It is assumed that the ribosome assembly step in protein translation, described by the k^ 
rate constant, is not rate limiting 



Observable value 


Initiation(fci) 


Step after 
initiation^) 


Ribosome assembly (ks) 


Elongation (fo^) 




Wild-type cap, inefficient initiation 




Steady-state rate 
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creases 
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no change 


decreases after 
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Relaxation time 


no change 


no change 


no change 


goes up and 
down 




Wild-type cap, efficient initiation 




Steady-state rate 


no change 


slightly de- 
creases 

after strong 
inhibition 


no change 
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Relaxation time 


no change 


goes up and 
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no change 


no change 




A 


-cap, inefficient initiation 




Steady-state rate 


decreases 


decreases 


no change 


slighly decreases 
after strong in- 
hibition 


Relaxation time 


no change 


no change 


no change 


goes up and 
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A-cap, efficient initiation 




Steady-state rate 


decreases after 
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slightly de- 
creases after 
strong inhibi- 
tion 


no change 


decreases 


Relaxation time 


goes up and 
down 


goes up and 
down 


no change 


increases 
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conditions are strongly different from situation in vivo. Indeed, inside the cell, mRNAs (microRNA 
targets) exist as ribonuclcoprotein particles or mRNPs, and second, all proteins normally associated 
with mRNAs transcribed in vivo are absent or at least much different from that bound to the same 
mRNA in an in vitro system or following the microRNAs transfection into cultured cells. The fact 
that RNA-binding proteins strongly influence the final outcome of microRNA regulation is proved now 
by several studies [31 IT2] 122] . 

The mathematical results provided in this paper suggests a complementary view on the co-existence 
of multiple microRNA-mediated mechanisms of translation repression. Mathematical modeling sug- 
gests to us to ask a question: if multiple mechanisms act at the same time, would all of them equally 
contribute to the final observable repression of protein synthesis or its dynamics? The dynamical limi- 
tation theory gives an answer: the effect of microRNA action will be observable and measurable in two 
cases: 1) if it affects the dominant system of the protein translationary machinery, or 2) if the effect 
of microRNA action is so strong that it changes the limiting place (the dominant system). 

In a limited sense, this means, in particular, that the protein synthesis steady rate is determined 
by the limiting step in the translation process and any effect of microRNA will be measurable only if 
it affects the limiting step in translation, as it was demonstrated in [37] . Due to the variety of external 
conditions, cellular contexts and experimental systems the limiting step in principle can be any in 
the sequence of events in protein translation, hence, this or that microRNA mechanism can become 
dominant in a concrete environment. However, when put on the language of equations, the previous 
statement already becomes non-trivial in the case of non-linear dynamical models of translation (and 
even linear reaction networks with non-trivial network structure). Our analysis demonstrates that the 
limiting step in translation can change with time, depends on the initial conditions and is not repre- 
sented by a single reaction rate constant but rather by some combination of several model parameters. 
Methodology of dynamical limitation theory that we had developed [16l [17] , allows to deal with these 
situations on a solid theoretical ground. 

Furthermore, in the dynamical limitation theory, we generalize the notion of the limiting step to 
the notion of dominant system, and this gives us a possibility to consider not only the steady state 
rate but also some dynamical features of the system under study. One of the simplest measurable 
dynamical feature is the protein synthesis relaxation time, i.e. the time needed for protein synthesis to 
achieve its steady state rate. The general idea of "relaxation spectrometry" goes back to the works of 
Manfred Eigen, a Nobel laureate [8] and is still underestimated in systems biology. Calculation of the 
relaxation time (or times) requires careful analysis of time scales in the dynamical system, which is 
greatly facilitated by the recipes proposed in [TTJ [35] . As we have demonstrated in our semi-analytical 
solutions, measuring the steady state rate and relaxation time at the same time allows to detect which 
step is possibly affected by the action of microRNA (resulting in effective slowing down of this step) . To 
our knowledge, this idea was never considered before in the studies of microRNA-dependent expression 
regulation. The table [2] recapitulates predictions allowing to discriminate a particular mechanism of 
microRNA action. 

Thus, analysis of the transient dynamics gives enough information to verify or reject a hypothesis 
about a particular molecular mechanism of microRNA action on protein translation. For multiscale 
systems only that action of microRNA is distinguishable which affects the parameters of dominant 
system (critical parameters) , or changes the dominant system itself. Dominant systems generalize and 
further develop the old and very popular idea of limiting step. Algorithms for identifying dominant 
systems in multiscale kinetic models are straightforward but not trivial and depend only on the order- 
ing of the model parameters but not on their concrete values. Asymptotic approach to kinetic models 
of biological networks suggest new directions of thinking on a biological problem, making the math- 
ematical model a useful tool accompanying biological reasoning and allowing to put in order diverse 
experimental observations. 
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